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Abstract 

This paper considers several topics arising in the finite element 
solution of the Incompressible Navier-Stokes equations. Specifically, the 
question of choosing finite element velocity/pressure spaces is addressed, 
particularly from the viewpoint of achieving stable discretizations leading to 
convergent pressure approximations. Following this, the role of artificial 
viscosity in viscous flow calculations is studied, emphasising recent work by 
several researchers for the anisotropic case. The last section treats the 
problem of solving the nonlinear systems of equations which arise from the 
discretization. Time marching methods and classical Iterative techniques, as 
well as some recent modifications are mentioned. 


Research reported here was suppported by the National Aeronautics and Space 
Administration under NASA Contract No. NASl-17130 while the authors were in 
residence at The Institute for Computer Applications in Science and 
Engineering, NASA Langley Research Center, Hampton, Va 23665. 
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Discretization 


Continuous Problen 

Let Q be a bounded region of I? or TP? , the flow region, and let 
u and p denote the velocity and pressure fields, respectively and v the 
kinematic viscosity- Normalizing the pressure by the constant density, the 
stationary Navier-Stokes equations take the form 

_u*Vu + Vp = vAu + ^ in (1) 

divu = 0 (2) 

_u = _g on dQ (3) 

where f and _g are given functions and dQ denotes the boundary of Q. 

In conservation form (1) may be written 

Div(u _u^ ) + Vp = vAu + f (4) 

where Div denotes the tensor divergence operator and the equivalence of (1) 
and (4) is shown using (2). Equation (2)-(4) are the equations to be solved. 
The standard weak form of (2)-(4) is: find satisfying (3) and 

2 

PGLq(Q) such that 

v/VurVv - /(n u*^] : Vv 

(5) 

- /pdlvv = /f*v V- veHQ(fi) 

/qdivu =0 V qeL^Cfi). 


( 6 ) 
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In the above formulation, H (f2) is the usual Sobolev space of functions with 

one square integrable derivative, is that subspace of whose 

2 

elements are zero on 9^2 and consists of square integrable functions 

with zero mean over Q. 

A slightly different treatment of the convection term is given in the 
following weak form of the momentum equation^ 

V / V^: V V - V 2 / (u • Vu* v-^* V v*_u) 

(7) 

- fpdivv = v£e1(Q). 

Q Q 

Again, in view of (2), this formulation is equivalent to (5); however for 
computations, the form (7) possesses certain advantages,^ and thus is 
recommended. 


Discrete Problem 

Proceeding in the usual way, one chooses finite dimensional subspaces 
V^CH^(fi), S^CLq(J 2), and seeks satisfying 

_u^ = on 9Q (8) 

and p^eS^ such that 

h „ h 1/ rr h „ h h h „ h h>i 

V J : Vv " V2 J [u .V^ 'Z. ” ^ J 

- Jp^dlvv^ = 

Q “ Q “ 


( 9 ) 
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/q^divu^ =0 ^ q^eS^. (10) 

In (8) denotes an approximation, in restricted to 9^2, to _g 

h h 

(e.g., an interpolant) . is the subspace of V of trial functions which 

are zero on Strictly speaking, this formulation is valid for polygonal 

domains but may be easily extended, e.g*, by isoparametric techniques, to more 
general domains. 

2 3 

Unlike the standard positive definite elliptic case * mere inclusion of 
^ ^ 2 

V CH (^2) and S CLq(? 2) is not sufficient to ensure convergence (or even 

Vi Vi 

existence) of the discrete solutions. In fact, the spaces and S” 

cannot be chosen independently of each other. Mathematically, the following 
condition is required^ 


sup 

V eV 
. h. 


h 

0 


V 


1 


/q^divv^ > yllq^ll 




( 11 ) 


where y > 0 is independent of the discretization parameter h. Here, the 
norms used in (11) are defined by 

Ivl^ = /Vv:Vv ^ 

llqll = Jq -V qeLQ(fi). 

n 

There are other mathematical conditions on the discrete spaces which must be 
imposed in order to guarantee convergence,^ however it is (11) which can fail 
to hold in general, and anyway offers substantial difficulties in its 
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verification. Equation (11) Is referred to as the condition of 
"dlv-stabllity".^’^ 

Note that the analogous condition 

sup /qdlvv > yllqll ^ qeLg(Sl) 

Ivl = 1 

is necessary to guarantee the existence and stability of the solution of the 
continuous problem, l.e., the Navier-Stokes equations. This condition Is 
easily verified since It Is well known that the problem 

dlvv * q In 

V = 0 on 9Q 

" — 1 2 

has a solution for any qzL^(Q) which satisfies 

yIvI < llqll 

for some constant y > 0. Then 

A A 

sup /qdlvv > /qdlvv / |v| 

- 

Ivl = 1 

= /q^ / Ixl ^ • 
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div-Stable Elements 

In practice, one would like to use low degree piecewise polynomial trial 
spaces V^, S^. Unfortunately these are the spaces which encounter the most 
difficulty in satisfying (11) • For the higher order cases, a reasonably 
efficient test for showing that an element pair is div-stable is known. ^ No 
such simple test is known for the low order cases. This is discussed further 
below. 

Perhaps the best known low order pair is the piecewise bilinear 
velocity /piecewise constant pressure combination defined on a quadrilateral 
subdivision of Q. This pair has been extensively used in engineering 
computations and been the object of much theoretical work. In particular, the 
checkerboard pressure mode is well documented.^ The effect of this mode is to 
make y in (11) zero, so (11) is caused to fail in this case. Less well 

o 

known is the existence of a mode such that y = 0(h). This causes the 
pressure approximation to fail to converge in general, even when the standard 
checkerboard mode is filtered out. This is proved in Boland and Nicolaides^, 
where also a filtering technique applicable to many cases is given. Use of 
this filter enables the pressures to converge optimally. 

Unfortunately, for arbitrary elemental subdivisions of Qy the 
appropriate filters are not known. Therefore, it seems advisable to use 
element pairs known a priori to be div-stable. Here are presented two element 
pairs which have been used in practical calculations and have been shown 
rigorously^ to be stable. 

For the first of these, let Q = { (x,y) ,0<x,y<l} and subdivide by 
lines X = ih, y = jh, i, j=0,l , . . . ,2n, so that h = l/2n. For V one may 
use continuous piecewise bilinear vector fields, bilinear in each subsquare. 
The pressure field is a subspace of all piecewise constants with zero mean on 
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defined as follows: on the coarser grid = 1 /n whose subsquares each 

contain four of the smaller subsquares, one merely Imposes one constraint on 
the four scalars. For example, referring to Fig. 1, one Imposes the 
constraint S]^ + S 2 = S 3 + s^ . For the general case. Isoparametric mapping 

onto this macroelement Is used. The resulting pair Is dlv-stable, requires no 
filtering and gives convergent approximations.^ 

Another stable pair, not requiring the use of Isoparametric methods Is 
the following. Let Q be triangulated regularly and let be an arbitrary 

element of the trlangulatlon. In t^, p^ is taken as a piecewise constant. 
To define Is further subdivided Into four similar triangles and 

Is then defined to be all continuous piecewise linear fields on the finer 
trlangulatlon. This element pair Is also dlv-stable and convergent without 
filtering. Fig. 2a shows a typical in the pressure trlangulatlon and 

Fig. 2b depicts the resulting velocity elements derived from x^. 

Since It Is known that these elements are dlv-stable, i.e., (11) holds, 

1 9 

then one may use finite element theory » to obtain the following estimates: 


Il_u-u^ll , < c-h 

E (Q) 


llp-p^ll ^ < c«h 

L (n) 


( 12 ) 


where c^^ and C 2 do not depend on h. The duality method then shows that 


llu-^^ll « < c«h^ . 

l"^(^) 


( 13 ) 


Thus, the methods are second order, in the root mean square sense, for the 
velocities and first order for the pressure. These are the best possible 
rates obtainable with these elements. 

Upwlndlng 

Since and are used as both test and trial spaces in (8)-(10) , 

the discrete equations will be of centered type. Thus, for fixed h as v^O 
one has to expect the numerical solution to develop increasingly large 
oscillations^^ ("wiggles”). These can of course be eliminated by use of any 
of the numerous upwinding/artif icial diffusion methods. However, the very 
idea of letting while keeping h fixed precludes the possibility of 

accurate computation of viscous effects. For this one must, at the very 
least, permit h to go to zero with v. The precise dependence of h on v 
depends on what viscous phenomena are to be accurately simulated. 

For example, simulation of shear or tangential boundary layers certainly 
will require h = ] in the neighborhood of these layers. Conceivably, 

especially at outflows, there may be other kinds of layers, maybe induced by 
numerical boundary conditions, requiring the more onerous (and less practical) 
restriction h =* 0(v). Usually the latter layers appear not to be of physical 
interest and may therefore be smoothed by use of anisotropic artificial 
diffusion, chosen so that tangential layers are not smoothed. Such methods 
have recently come into prominence. Specifically these methods attempt 
to add diffusion only in the streamwise direction. One must emphasize that 
h = o[v^^ ) is still an essential requirement for resolution of physical 
layers and hence, convergence. 
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One streamwlse artificial diffusion method is now presented. First, 
observe that (4) may be written in the form 

Dlv[_u u^+pI-vV_u] = (14) 

where I is the Identity matrix. The quantity in the brackets is the 

momentum flux density. To this one adds the term 

- AVu(u (15) 

where X is a mesh dependent parameter tending to zero with h. Then (14) 
becomes 

Dlv[u u^+pI-V^[ vI+X^ ] = _f (16) 

In the form (16) one may Interpret (15) as an anisotropic perturbation of the 

viscosity tensor (vl). Clearly the effect of this perturbation vanishes in 

directions normal to the velocity vector, and hence to streamlines. 

Alternately, by associating the perturbation with the convection term a 

connection may be made, in transient cases, with the method of 

characteristics . This is not pursued further here . Equation ( 16) , along 

with (2)-(3), may now be discretized in the usual manner using the same test 

and trial spaces. Generally, in order to achieve the desired stabilizing 

effect, X should be 0(h). The effect of such a choice of X on the 

estimates (12)-(13) is not fully understood. When using higher order elements 

it is reasonable to expect that the accuracy would be degraded by this crude 

approach. Therefore it is of Interest to note that a similar effect may be 

1 3 

achieved by approplate choices of distinct test and trial spaces. 
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Normally, for internal flows, e.g. , cavity flows, and for good choices 
for numerical outflow conditions, 0(v) layers are not present in the flow 
field. In such cases, the use of artificial streamwise diffusion methods is 
not necessary. 

One is still left with the task of resolving tangential layers. If a 
uniform mesh is used throughout the flow field, the h = o(v^2 J restriction 
results in an unacceptable number of degrees of freedom. However such a 

small h is needed only in the neighborhood of the layers themselves. Thus 
it is advantageous, in the sense of reducing the number of degrees of freedom, 
to use nonuniform grids. As an example consider the driven cavity problem. 
Here Q is a unit square, the upper end of the cavity moves with velocity 
- (1,0), and V = 1/3200. The results of three calculations are reported. 
The first^^ uses an upwind finite difference technique, on a uniform mesh, for 
a streamfunctlon - vorticity formulation of (1) - (3). The second^^ again 
solves the streamfunctlon - vorticity formulation by a finite element 
technique using a nonuniform grid. The grid spacing is determined by the 
functions 

2 — 

X = sin (ttx/2) 

(17) 

2 _ 

y = sin (7Ty/2) 

with a uniform grid spacing in the x and y coordinates. The third 
calculation uses the element pair of Fig. 2 in conjunction with the primitive 
variable formulation and a nonuniform grid again determined by (17). In 
neither finite element calculation were artificial dif fusion/upwlndlng 
techniques used. In Fig. 3 is given the y-component of velocity v at 
y = 1/2 as a function of x. The uniform grid calculation used a 
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(129 X 129) grid so that h = 0.0078. The nonuniforra grid calculations used 
a (19 X 19) grid. With the spacing determined by (17), the minimum 
h = 0.0076 which is comparable to the uniform h of the calculation of Ghla, 
et.al..^^ Clearly, one can achieve the same accuracy at a greatly reduced 
cost by using nonuniform grids. 


Solution Technique 

Whatever method is employed for discretization, the outcome is a large 
system of nonlinear equations which must be solved for the approximate 
solution. Concerning the solution of this system by some classical technique, 
such as Newton^s method, one must first observe that the Jacobian of the 
system requires a remarkably large amount of storage, if the usual band 
storage scheme is adopted. For example, it is easy to verify that a driven 
cavity calculation on an n x n grid requires 27n^ words of storage (for the 
bilinear/constant element). Thus, even with n = 20, a rather large storage 
requirement is apparent. In 3D, the situation is catastrophic. Hence, the 
classical techniques may be of limited value (although some fairly fine mesh 
calculations have been reported^^). Anyway, it is evident that alternative 
techniques are of interest. Of course many have been proposed. They fall 
into the two natural classes of "false transient" or time marching algorithms 
on the one hand and general nonlinear equation solvers on the other. 

Here is an example of a time marching technique used successfully by the 
authors. Let U, P denote the unknowns and U^,P^ the k-th approximations 
to U and P, and consider 
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(18) 


k - 0,1,*** 


BTyk+l 


= G 


(19) 


k 


0 , 1 , 


• • • 


where N(U) represents an approximate Laplacian plus convection terms and 
B an approximate gradient operator, generated by the particular finite 
elements used. The first term In (19) Is a discretization of -|v, forward In 

O L 

time from the time level. is arbitrary, and is not required to 

satisfy (18) . In the steady state limit, and one recovers the 

solution of the original problem. To get , given U^, multiply (18) by 

and use (19) to get 


BlBpk+l ^ gTp _ Vi > I 


k+1 

which must be solved for BP , so that then 


If k = 0, a term 


yk+l ^ yk At[F-BP^‘^^-N(u’^) ] 




( 20 ) 


( 21 ) 


is added to the right hand side of (20) • Notice that (20) is the usual 
"Poisson" equation for the pressure, but with the crucial feature that 
boundary conditons do not need to be supplied externally. They are 
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automatically present in . Notice also that in (21) not 

k+1 

BP^ is required • This is the least squares solution of the equation 


where 


gTvk+l 


b^’f - b^n[u^) 


,k+l 


BP 


,k+l 


( 22 ) 


which may be found by any of numerous techniques, direct or iterative. The 
following Iterative method, a generallztlon fo a method of Kacmarcz which may 
be called the row projections method, was used by the authors. For solving 


T 

B V = L 


this method is the following: 


£+1 £ 
VXTl = yX, ^ 


V" = B?, 


arbitrary, £=1 ,2 , • 


where B, 


denotes the 


,th 


column of B (taken cyclically), o) is a 

,th 


relaxation parameter and a a number chosen so that when o) = 1, the Z 

£+1 T £+1 

component of the residual r = L - B V is zero, thus, 


a = 


2 T 

where is the sum of the squares of the elements in row £ of B . 

Geometrically, the error is projected orthogonally in succession 

onto the planes whose noxrmals are the rows of B^. o) effects a kind of 
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overprojection. Generalization to projection onto the planes normal to 
several rows (corresponding to line, etc., relaxation) is straightforward. 

Returning to (20), note that only few iterations of (22) are required for 
each time step, due to the availability of the solution from the earlier time 
level. 

The problem with time marching based on first order time derivatives is 

the large number of time steps (regardless of the size of v) needed to reach 

the steady state. Generally, several thousand steps are needed for each digit 

of accuracy required, for the kind of grid sizes encounted in practice, and 

the number of steps rapidly Increases as the latter approach zero. Using 

higher order time derivatives would give considerable improvement, but does 

not appear to be a common idea. Among the other methods, great promise is 

17 18 

shown by the reduced basis - continuation techniques . ^ A nonlinear 

1 q 

conjugate gradient method is used by Glowinski, et.al. , although it does not 
appear to function well in all cases. Various iterative techniques have also 
been applied with mixed success by the authors . These will be reported 
elsewhere. 


Conclusions 

In connection with the incompressible Navier-Stokes equations 
discretization techniques via finite elements, accuracy questions, and methods 
for solving the algebraic systems of nonlinear equations, all for the u ,p 
(primitive) variable case, have been discussed. It is fair to say that the 
first topic is by now reasonably well understood. The main difficulty 


concerned the discretization of the incompressibility condition, and this 
topic now has a variety of theoretical analyses and tests whereby the 
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stability of the discretization can be verified before calculations are 
performed. 

Less satisfactory is the state of affairs relating to solving the systems 
of nonlinear equations. Here, there are many problems still to be overcome, 
mostly concerning the efficiency of the solution procedures. The time 
marching method discussed, along with similar methods based on first order 
time derivatives really is too inefficient, requiring of the order of lO*^ - 
10^ time steps to obtain each digit of the solution. On the other hand, the 
direct solution methods encounter problems caused by the large storage 
resources necessary for carrying out the matrix manipulations as well as 
domain of convergence problems as v ->• 0. In the latter case, the starting 
approximation has to be closer and closer to the exact solution being 
computed, as v becomes smaller. Continuation techniques are naturally 
suggested for dealing with this latter issue, as stated in the text. 
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